home *** CD-ROM | disk | FTP | other *** search
/ Almathera Ten Pack 3: CDPD 3 / Almathera Ten on Ten - Disc 3: CDPD3.iso / scope / 201-220 / scopedisk202 / bbbbs2 / rexxdoors / moon.rexx < prev    next >
OS/2 REXX Batch file  |  1995-03-19  |  4KB  |  156 lines
  1. /* Moon.rexx        by Greg Cunningham 1990       converted from Unix TODAY.C
  2.  *
  3.  * GMT Settings: (America)
  4.  *
  5.  * Eastern  5
  6.  * Central  6
  7.  * Mountain 7
  8.  * Pacific  8
  9.  *
  10.  */
  11.  
  12. CR='0D'x
  13. TimeZone=8
  14. DaylightSavings=0 /* set to 1 if daylight */
  15.  
  16. if ~Show(L,'rexxmathlib.library') then
  17.    if ~AddLib('rexxmathlib.library',0,-30,0) then
  18.       Say 'No RexxMathLib, running anyway.'CR
  19.  
  20. say lunar(date('S'),time(),TimeZone-DaylightSavings)||CR
  21.  
  22. exit
  23.  
  24. /* Phase of the Moon. Calculates the current phase of the moon.
  25.  * Based on routines from `Practical Astronomy with Your Calculator',
  26.  * by Duffett-Smith.
  27.  *
  28.  * Synopsis:
  29.  *
  30.  * string=lunar(year,time,gmt)
  31.  * date  -- Current sorted date
  32.  * time    -- Current time of the day
  33.  * gmt     -- Hours from GMT time zone (subtract one if daylight savings)
  34.  *
  35.  * say lunar(date('S'),time(),0) current phase of the moon
  36.  * say lunar(date('S')+1,time(),0)  same time tommorow
  37.  *
  38.  * say lunar('19900225','12:00:00',0)  Feb 25,1990 at noon
  39.  *
  40.  */
  41.  
  42. lunar:   procedure
  43.    arg cdate,ctime,GMToffset
  44.  
  45.    /* days between Jan 1,1985 and Now */
  46.    days=date('B',cdate,'S')-date('B','19850101','S')+1
  47.  
  48.    /* add fraction of day */
  49.    parse var ctime hour':'min':'sec
  50.    hour=hour+GMToffset
  51.    if(hour>23) then do
  52.       days=days+1
  53.       hour=hour-24
  54.    end
  55.    /* add on percentage of day */
  56.    days=days+((hour+(min/60)+(sec/3600))/24)
  57.  
  58.    phase=potm(days)
  59.  
  60.    cp='The Moon is'
  61.  
  62.    if(phase=0) then     cp=cp 'New'
  63.    else if(phase=100) then    cp=cp 'Full'
  64.    else if(phase=50) then do
  65.       phase2=potm(days+1)
  66.       if(phase2>phase) then   cp=cp 'at the First Quarter'
  67.       else        cp=cp 'at the Last Quarter'
  68.    end
  69.    else if(phase>50) then do
  70.       phase2=potm(days+1)
  71.       if(phase2>phase) then   cp=cp 'Waxing'
  72.       else        cp=cp 'Waning'
  73.       cp=cp 'Gibbous ('phase'% Full)'
  74.    end
  75.    else do
  76.       phase2=potm(days+1)
  77.       if(phase2>phase) then   cp=cp 'Waxing'
  78.       else        cp=cp 'Waning'
  79.       cp=cp 'Crescent ('phase'% Full)'
  80.    end
  81.    return cp
  82.  
  83.  
  84. /* calc phase-of-the-moon; returns percentage of full */
  85. potm: procedure
  86.    arg days
  87.  
  88.    EPSILONg=279.611371  /* solar ecliptic long at EPOCH */
  89.    RHOg  =282.680403 /* solar ecliptic long of perigee at EPOCH */
  90.    e  =0.01671542 /* solar orbit eccentricity */
  91.    lzero =18.251907  /* lunar mean long at EPOCH */
  92.    Pzero =192.917585 /* lunar mean long of perigee at EPOCH */
  93.    Nzero =55.204723  /* lunar mean long of node at EPOCH */
  94.    PI =3.141592654
  95.  
  96.    N= 360 * days / 365.2422   /* days*360/real_days_of_a_year */
  97.    N=pos360(N)       /* round up to positive 360 degrees */
  98.  
  99.    Msol=N+EPSILONg-RHOg
  100.    Msol=pos360(Msol)
  101.  
  102.    Ec=360 / PI * e * sin( dtor(Msol) )
  103.  
  104.    LambdaSol = N + Ec + EPSILONg
  105.    LambdaSol = pos360(LambdaSol)
  106.  
  107.    l = 13.1763966 * days + lzero
  108.    l = pos360(l)
  109.  
  110.    Mm = l - (0.1114041 * days) - Pzero
  111.    Mm = pos360(Mm)
  112.  
  113.    Nm = Nzero - (0.0529539 * days)
  114.    Nm = pos360(Nm)
  115.  
  116.    Ev = 1.2739 * sin( dtor(2*(l - LambdaSol) - Mm) )
  117.  
  118.    Ac = 0.1858 * sin( dtor(Msol) )
  119.    A3 = 0.37 * sin( dtor(Msol) )
  120.  
  121.    Mmprime = Mm + Ev - Ac - A3
  122.  
  123.    Ec = 6.2886 * sin( dtor(Mmprime) )
  124.  
  125.    A4 = 0.214 * sin( dtor(2 * Mmprime) )
  126.  
  127.    lprime = l + Ev + Ec - Ac + A4
  128.  
  129.    V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol)))
  130.  
  131.    ldprime = lprime + V
  132.  
  133.    D = ldprime - LambdaSol
  134.  
  135.    perc= 50 * (1 - cos( dtor(D) ))
  136.  
  137.    return trunc(perc+.5)   /* return percent rounded to nearest integer */
  138.  
  139.  
  140.  
  141. /* get a positive 360 degree mod */
  142. pos360: arg deg
  143.  
  144.    deg=deg//360
  145.    if(deg<0) then deg=deg+360
  146.    return deg
  147.  
  148.  
  149. /* convert degrees to radians */
  150. dtor: procedure
  151.    arg deg
  152.  
  153.    PI=3.141592654
  154.    return (deg * PI) / 180
  155.  
  156.